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Abstract 

We study zero-temperature properties of the 3d Edwards-Anderson Ising spin glass 
on finite lattices up to size 12 3 . Using multicanonical sampling we generate large num- 
bers of groundstate configurations in thermal equilibrium. Finite size scaling with a zero- 
temperature scaling exponent y = 0.74 ± 0.12 describes the data well. Alternatively, a 
descriptions in terms of Parisi mean field behaviour is still possible. The two scenarios give 
significantly different predictions on lattices of size > 12 3 . 
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The theoretical problem to determine the equilibrium groundstate structure of spin 
glasses has remained an important, but elusive, question. It is generally agreed that the 
statistical mechanics of the infinitely ranged Sherrington-Kirkpatrick Ising spin glass is es- 
sentially understood. The replica-symmetry breaking mean field (MF) scheme discovered 
by Parisi [jlj exhibits infinitely many low-temperature states whose properties are consistent 
with simulations 0, see || for a recent overview. Whether, for temperatures below the 
spin glass phase transition point T c , the more realistic short-ranged 3d Edwards-Anderson 
Ising spin glass model (EAI) also exhibits Parisi mean field behaviour has become a central 
question. It was answered in the negative by proponents of a simple scaling ansatz [4-7]. 
These droplet scaling (DS) theories suggest that no more than two pure states (related via 
a global flip) exist at any temperature. The MF approximation is surely valid for d — ► oo, 
and it has been suggested H that d — 6 is the upper critical dimension which separates the 
MF from the DS scenario. 

The EAI Hamiltonian is given by 

H = - J ii s i s i- (!) 

<ij> 

Here the sum < ij > goes over nearest neighbours. We consider 3d systems with periodic 
boundary conditions and N = L 3 spins. The exchange interactions Jy = ±1 between N 
spins Si = ±1 are randomly distributed over the lattice with the constraint J2<ij> J%j — 0. 
For each system there are (3A)!/[(3A/2)!] 2 realizations of the quenched random variables 
J = {Jij}- Recent simulations || in a magnetic field favour the mean field picture rather 
than the alternative droplet model. However, it has been pointed out that equilibrium at 
sufficiently low temperatures has not been reached [[|. 

A quantity of decisive importance is the probability density P(q) of the Parisi order 
parameter q: 

p (?H< p W>J .ra£P j(5) . (2) 

By < ■ >j we denote averaging over the realizations J. For a fixed realization Pj(q) is the 
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probability density of the overlap 
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Here and s 2 denote two replica (i.e statistically independent configurations) of the realiza- 
tion J at temperature T. Due to magnetic field zero we have the symmetry P{—q) = P{q), 
such that J^l q n P(q)dq = for n odd. We therefore define averages over the range < q < 1: 

=< g n >= 2 f 1 q n P(q)dq. (4) 
JO 

For the spin glass susceptibility Xq = N < q 2 > the MF as well as the DS scenario suggest 
divergence Xq ~ N for T < T c . However they give significantly different predictions for the 
variance 

a\q)=<{q-qf> . (5) 

In the limit iV — > oo one has er 2 (g) — * finite in MF theory, while cr 2 (g) ~ L~ y — > within 
the DS approach. Here ?/ = — j/t is the zero-temperature scaling exponent 0, denoted # in 
JJ, which governs also the finite size scaling (FSS) corrections of the expectation values 1 : 
g n L — q 71 ^ ~ for DS, whereas we assume 1/Volume corrections for MF theory. 

Lack of self-averaging is one prominent feature of MF behaviour. In reference || 
< f(P(q) — Pj(q)) 2 dq >j was studied. This was criticised by the authors of 0. Following 
their suggestion we estimate 

aj(q^)=<(7-7j) 2 >J- (6) 

Again, this quantity stays finite in MF theory, but drops off ~ L~ y in the DS picture. 

With the development of multicanonical techniques for disordered systems |10| it has 



become feasible to generate spin glass groundstates in thermal equilibrium; see [|ll|] for a 



brief, general review, and |12| for the earlier umbrella sampling. A pilot study for the 



model at hand has been presented in [ 13[ . In essence a multicanonical spin glass simulation 



proceeds in three steps. First Monte Carlo (MC) weights are recursively constructed which 
1 If necessary we indicate the lattice size by an additional subscript L, which is otherwise dropped. 



will allow to simulate an ensemble, the "multicanonical" , which yields canonical expectation 
values in the temperature range < T < oo through use of the spectral density. Secondly, 
equilibrium configurations with respect to the multicanonical ensemble are generated by 
means of standard MC. In a third step canonical expectation values at desired temperatures 
are obtained from the analysis. Multicanonical sampling circumvents the notorious ergodicity 
problems of canonical low temperature spin glass simulations through regular excursions into 
the disordered phase, while staying in equilibrium. 

In this paper we focus on the investigation of groundstate properties. A lower bound 
on the number of statistically independent groundstates sampled is obtained by counting 
how often the system moves from the energy E > region to the groundstate energy E m i n , 
and back to the E > region. This has been termed "tunneling" |T0| and we follow this 
notation, but one should bear in mind that the free energy barriers are actually not overcome 
by a tunneling process. With present techniques the tunneling time approximately increases 
~ V in 3d [|T^]. This slowing down limits our investigation to rather moderately sized 
lattices. 

We have performed simulations for L = 4, 6, 8, 12 (N = 64, 216, 512, 1728). For L < 8 
the sum (2) is approximated through 512 randomly chosen realizations of the {Jij}, whereas 
we have only 7 realizations for L = 12. For all 1,543 cases multicanonical parameters were 
determined recursively. Then each system was simulated twice with independent random 
starts and random numbers. This constitutes our two independent replica per realization. 
In these production runs iterations were stopped when a preset number of tunneling events 
n T had occurred: n T = 128 (L = 4), 64 (L = 6), 32 (L = 8) and 10 (L = 12). Despite 
this decrease in tunneling events, the average number of updates per spin n s (sweeps) did 
steadily increase. Approximate values are: n s = 8- 10 4 (L = 4), 10 5 (L = 6), 7.6 • 10 6 (L = 8) 
and 50 • 10 6 (L = 12). The average CPU time spent on one L = 8 replica was approximately 
800 minutes on an IBM 320H workstation. 

Per replica we have stored up to 2,048 groundstate configurations. Due to correlations 
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the number of encountered groundstates is, of course, much larger than n T . If the number 
exceeded 2,048, the stored configurations were randomly selected from the total set. For 
groundstate configuration n this is elegantly done on-line by picking a random integer i r 
in the range 1 < i r < n. Configuration n is stored at position i r if i r < 2, 048 and dis- 
carded otherwise. For both replica the same groundstate energy has to be found between all 
tunneling counts. This is a strong, albeit not rigorous, criterium to ensure that the correct 
groundstate energy has not been missed. 

On a semi-log scale figure 1 depicts the thus obtained probability densities (2) for the 
Parisi order parameter. The L = 12 probability density, presented without error bars, is 
very bumpy due to the small number of realizations, and will only be reliable for a few 
of the considered physical quantities. Note, altogether the data respect the P(q) = P{—q) 
symmetry well. For L = 8 figure 2 plots the Pj{q) probability densities of two rather extreme 
L = 8 realizations: two peak shape versus continuous distribution. Various different shapes 
in-between these extremes are also found. Figure 3 shows all L = 8 realizations together. 
From figures 1-3 it is evident that only a careful quantitative analysis of these distributions 
may give hints concerning the L — > oo groundstate distribution. 

Our estimates for various measured quantities are summarized in table 1. The error 
bars are with respect to the different realizations, which are statistically independent and 
enter with equal weights. Table 2 summarizes two-parameter fits of the data, assuming 
alternatively MF theory or the DS ansatz to be true. If MF and DS scenario lead to the 
same functional form, the fit is marked "All". The oo column gives the infinite volume 
extrapolations of the considered quantity, Q is the goodness of fit, and R 12 comments Yes 
or No on the reliability of the L — 12 data for the purposes of the particular fit. With an 
exception for P max , the MF and ALL fits are of the form a\ + a 2 /L 3 . In the MF as well as 
in the DS scenario the self-overlap gives rise to a ^-function singularity. Therefore P max is 
supposed to grow ~ V and the appropriate fit is aiL 3 + a 2 . Including all data points the fit 
is still consistent. Although, omitting the smallest lattice indicates that L = 4 may not fully 
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exhibit the asymptotic behaviour. 

Let us now discuss a 2 (q) and cr 2 (q 2 )- The DS fit is a\L~ a2 . Fits and data are depicted 
in figure 4. For each case we give two fit curves. The upper one relies on three data points 
(L = 4,6,8), whereas the lower one includes also the L = 12 result. When only three data 
are used, MF and DS fits are both consistent (Q = 0.10 and Q = 0.43). Once the L — 12 
data point is included, the consistency of the MF fit becomes marginal (Q = 0.04). However, 
from the L = 8 data we have the experience that 10% of the realizations amount to 99% of 
the P(0) contribution. Consequently, the L = 12 data suffer not only from large statistical 
fluctuations, but are altogether unreliable for quantities which are sensitive to the small q 
distribution. The above fits were also used for -P(O), but the data are too inaccurate to yield 
meaningful results. We now rely on the three-point fits for <r 2 (q) and a 2 (q 2 ). The two y = a 2 
estimates are still compatible, and we summarize them to y = 0.74 ± 0.12. The error bar is 
not reduced as both estimates rely on the same data set. 

Assuming that our L = 4 — 8 lattices show already typical scaling behaviour, we conclude 
from figure 4 that similarly accurate data on lattices up to size L — 16 would discriminate 
between the MF and the DS ansatz. Due to the slowing down, our new data indicate ~ V 3S , 
the needed CPU time would be about 1,000 times larger than the one spent on the present 
investigation. With upcoming massively parallel devices in the teraflop range such a factor 
can be achieved. 

The not yet discussed DS fits are of the form a\ + a 2 L~ y with y = 0.74. Let us first 
comment on the groundstate energy fits. The energy is self-averaging, L — 12 contributes 
reasonably accurate results, and we rely on all data sets. Most interesting, the DS fit is 
consistent, whereas the MF fit with 1/Volume corrections is ruled out. Unfortunately, there 
is still a catch to it. It may well be that the corrections to the uncritical MF behaviour 
are exponentially small. Although FSS corrections for larger systems would then be greatly 
reduced, the disadvantage at the present level is that the appropriate a± + a 2 exp(— a 3 L) fit 
has three free parameters. A consistent (Q = 0.19) fit is then still possible. Again, accurate 
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data on lattices up to L = 16 would allow to differentiate this behaviour from DS. 

Let us remark that the groundstate entropy stays finite, as about 5% of the spins can 
still be freely flipped due to the exact degeneracy of the ±1 quenched random variables [14 |. 



Because of this pathology of the model non-critical FSS correction seem to be appropriate 
in either scenario, and 1/Volume corrections work indeed well. 

A relevant consistency check for the correctness of the DS picture is that the infinite 
volume estimates of <? max , \J~~cf and q have to agree. Figure 5 shows the MF and DS fits for 
these quantities. For L a log-scale is used to exhibit L —>■ oo clearly, and the infinite volume 
estimates are depicted towards the end off the scale. With the previously determined zero- 
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temperature exponent the values are indeed consistent. For the MF fits q^ < y q 2 ^ < q 
as it should be then. It is notable that fitting with a wrong zero-temperature exponent may 
produce entirely inconsistent results. For instance with y = 0.2 one finds q > 1 > q^ax- 

In summary, the DS ansatz is so far consistent. Our investigation presents the first MC 
estimate of the zero-temperature scaling exponent. Obviously, our lattices are too small to 
allow seminal results. In particular, the MF picture is still a valid alternative. It is clear 
that, either by brute computer power or by algorithmic improvements, simulations on larger 
lattices will become feasible. It seems, we are approaching a numerical conclusion about the 
correct groundstate picture of the 3d EAI model. 

Acknowledgements: We would like to thank David Huse for a useful e-mail communi- 
cation. Our simulations were performed on the SCRI cluster of fast workstations. 
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Table 1: Data. 
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Figure Captions 



Fig. 
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Probability densities P(q) for the Parisi order parameter (L = 4,6,8 and 12) 


Fig. 


2 


Probability densities Pj{q) for two very different L = 8 realizations. 


Fig. 


3 


All Pj(q) probability densities for L — 8. 


Fig. 
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Fits for (T 2 (q) (Data q) and crj(q 2 ) (Data q2). 


Fig. 
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Fits and extrapolations for g max (up), \f^f (middle) and q (down). 
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Table 2: Fits. 
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